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ABSTRACT 

<n : 

\ We study the hydrodynamical evolution of massive accretion disks around 

black holes, formed when a neutron star is disrupted by a black hole in a binary 
system. The initial conditions are taken from three-dimensional calculations of 
coalescing binaries. By assuming azimuthal symmetry we are able to follow the 
time dependence of the disk structure for 0.2 seconds in cylindrical coordinates 
• \ (r, z). We use an ideal gas equation of state, and assume that all the dissipated 

energy is radiated away. The disks evolve due to viscous stresses, modeled with 
an a-law. We study the disk structure, and in particular the strong meridional 
circulations that are established and persist throughout our calculations. These 
consist of strong outflows along the equatorial plane that reverse direction close 
to the surface of the disk and converge on the accretor. In the context of gamma 
ray bursts (GRBs), we estimate the energy released from the system in neutrinos 
and through magnetic-dominated mechanisms, and find it can be as high as 
E v ph 10 52 erg and Ebz ~ 10 51 erg respectively, during an estimated accretion 
timescale of 0.1-0.2 seconds, vv annihilation is likely to produce bursts from only 
a short, impulsive energy input L vV oc t -5 / 2 and so would be unable to account 
for a large fraction of bursts which show complicated light curves. On the other 
hand, a gas mass ~ 0.1 — O.25M survives in the orbiting debris, which enables 
strong magnetic fields ~ 10 16 G to be anchored in the dense matter long enough 
to power short duration GRBs. We highlight the effects that the initial disk and 
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black holes masses, viscosity and binary mass ratio have on the evolution of the 
disk structure. Finally, we investigate the continuous energy injection that arises 
as the black hole slowly swallows the rest of the disk and discuss its consequences 
on the GRB afterglow emission. 

Subject headings: accretion, accretion disks — hydrodynamics — gamma rays: 
bursts 



1. Introduction 

Accretion onto black holes has been considered as an efficient way to transform gravita- 
tional energy into radiation (Salpeter 1964; Zel'dovich 1964), and is often thought to occur 
in the form of a disk, due to the angular momentum of the accreting matter. This almost 
certainly is the case in a variety of astrophysical systems, ranging from AGNs with very mas- 
sive black holes to stellar mass binaries, analogous to the X-ray binaries known to contain 
neutron stars (King 1995). The flows in these disks may exhibit very different morphologies 
depending on the physical conditions present in the system. Accretion is generally thought 
to proceed through the transport of angular momentum from the inner to the outer regions 
of the disk, although the mechanism by which this is accomplished is not entirely clear. The 
parametrization introduced by Shakura & Sunyaev (1973) has allowed much progress to be 
made, without specifying the physics behind the viscosity responsible for angular momentum 
transport. Magnetohydrodynamical studies (analytical and numerical) appear to indicate 
that magnetic fields and their associated instabilities in disks can effectively generate a vis- 
cosity that would drive their evolution (Balbus & Hawley 1991; Hawley & Balbus 1991; 
Hawley 2000; Stone & Pringle 2001; Hawley & Krolik 2001, 2002), with equivalent values 
for the a parameter in the range 0.01-0.1 (Balbus & Hawley 1998). Additionally, it has 
become clear that hydrodynamic processes can play an important role in the structure and 
evolution of disks, and that multi-dimensional, time-dependent computations are necessary 
to fully understand these effects, particularly since the flows can be quite complicated, often 
exhibiting strong variability, and composed of combinations of inflows/outflows of varying 
intensity (Igumenshchev 2000). 

In this context we are motivated to study accretion disks around black holes hydrody- 
namically, particularly in what concerns central engines for cosmological gamma ray bursts 
(GRBs). The energetics of GRBs (10 50 — 10 52 erg are typically released in a few seconds) and 
the variability shown in their lightcurves (down to millisecond timescales) argues in favor of 
a compact source that produces a relativistic outflow, usually referred to as a fireball (Rees 
& Meszaros 1992). The complicated light curves can then be understood in terms of internal 



- 3- 



shocks in the outflow itself, caused by velocity variations in the expanding plasma (Rees & 
Meszaros 1994; Sari & Piran 1997; Fenimore, Ramirez-Ruiz & Wu 1999; Ramirez-Ruiz & 
Fenimore 2000; Nakar & Piran 2002). The durations range from 10~ 3 s to about 10 3 s, with 
a bimodal distribution of short (~ 0.2 s) and long bursts (~ 40 s) (Kouveliotou et al. 1993; 
Norris et al. 1996). One possible scenario, at least for the class of short bursts 1 , involves 
the coalescence of compact binaries containing a black hole (BH) and a neutron star (NS), 
or two neutron stars (Lattimer & Schramm 1976; Paczyhski 1986, 1991; Eichler et al. 1989; 
Narayan, Paczyhski & Piran 1992; Mochkovitch et al. 1993; Katz & Canel 1996; Kluzniak & 
Lee 1998; Popham, Woosley & Fryer 1999; Ruffert & Janka 1999; Rosswog & Davies 2002). 
Such systems do exist, like PSR1913+16 (Hulse & Taylor 1975), and they will coalesce due to 
angular momentum losses to gravitational radiation, provided that the orbital separation is 
small enough. The binary coalescence rates are in rough agreement with the observed GRB 
rate (Kalogera et al. 2001). The typical dynamical timescale in such binaries immediately 
prior to coalescence (~ ms) is much shorter than the observed burst duration, and so it 
requires that the central engine evolves into a configuration that is stable, while retaining a 
sufficient amount of energy to power the burst. 

The formation of a BH with a debris torus around it is a common ingredient of both these 
scenarios, whose accretion can provide the release of sufficient gravitational energy ps 10 54 
erg to power a GRB (Rees 1999). A fireball arises from the large compressional heating and 
dissipation associated with this accretion, which can provide the driving stress necessary for 
relativistic expansion (Piran 1999; Meszaros 2001). Possible forms of this outflow are kinetic 
energy of relativistic particles generated by vv annihilation or an electromagnetic Poynting 
flux. 

In either mechanism, the duration of the burst is determined by the viscous timescale 
of the accreting gas 2 , which is significantly longer than the dynamical timescale, thus ac- 
counting naturally for the large difference between the durations of bursts and their fast 
variability. Any instability (of hydrodynamical or magnetic origin) would presumably be 
reflected in the relativistic outflow. Strong magnetic fields anchored in the dense matter 
surrounding the BH would produce large amplitude variations in the energy release (Usov 
1992; Meszaros & Rees 1997). A weaker field would extract inadequate power; on the other 
hand a neutron torus, with its huge amount of differential rotation, is a natural site for 
the onset of a dynamo process that winds up the magnetic field to the required intensity 



1 We note, furthermore, that those with detected afterglows are all in the long category (although see 
Lazzati, Ramirez-Ruiz & Ghisellini (2001) and Connaughton (2002)). 

2 In the collapsar scenario, the burst duration is given by the fall-back time of the gas (Woosley 1993; 
MacFadyen & Woosley 1999) 
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(Usov 1994; Kluzniak & Ruderman 1998). An acceptable model requires that the orbiting 
debris not be dispersed completely on too short a timescale, lasting at least as long as the 
characteristic duration of the burst. The tidal disruption of a neutron star by a black hole 
would redistribute matter and angular momentum very rapidly. The key issue is then how 
long a sufficient amount of this matter survives to power a burst. 

In this paper, we study the evolution of realistic disks resulting from dynamical coales- 
cence calculations (see below §2), on timescales that are comparable to the durations of short 
GRBs (i.e. a few tenths of a second). Recently, the steady state structure of similar disks has 
been examined by Narayan, Piran & Kumar (2001) and Kohri & Mineshige (2002). Since we 
wish to investigate the evolution of the disk and its stability, no steady state assumptions are 
made regarding its structure. This allows us to compare the strength of the energy release 
by vv annihilation relative to MHD coupling and to investigate their potential as a viable 
source for GRB production. 

The hydrodynamics of black hole-neutron star coalescence has been considered in detail 
by Lee & Kluzniak (1999a,b) and Lee (2000, 2001a), using an ideal gas equation of state and 
varying its stiffness through the adiabatic index T. The simulations explored the effect of 
different initial binary mass ratios q b = M^s/Mbr and the spin configuration of the neutron 
star, with respect to the orbital motion. The viscosity inside neutron stars is far too small to 
permit synchronization during the inspiral phase (Kochanek 1992; Bildsten & Cutler 1992), 
so it is reasonable to assume that the neutron star spin is negligible with respect to the orbital 
angular momentum. The exact equation of state at supra-nuclear densities is uncertain, but 
it would appear that the radius is largely independent of the mass (Lattimer & Prakash 
2001), so that in the polytropic approximation, one would assume T = 2. Coalescence 
calculations with this set of parameters were computed by Lee (2001a, herafter L01). In 
particular, we will use here the results of runs C50 and C31 in that paper, with q b = 0.5 and 
q b = 0.31 respectively (both runs used an index T = 2). 



2. Numerical method and implementation 

The equations of hydrodynamics are solved in two dimensions using Smooth Particle 
Hydrodynamics (SPH) (Monaghan 1992). Three-dimensional dynamical coalescence calcu- 
lations are only followed for a few tens of milliseconds, at most (those in L01 spanned 23 ms), 
because of computational limitations. Since we wish to explore the behavior of the accretion 
structures on timescales that are much longer than the dynamical one, we have assumed 
azimuthal symmetry in the system, and mapped the output from the 3D calculations to 
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cylindrical (r, z) coordinates. This allows us to follow the calculations for a few tenths of a 
second (all the runs shown in this paper were terminated at t — 0.2 s). We have not assumed 
reflection symmetry with respect to the equatorial (z = 0) plane. 

We use a simple ideal gas equation of state, with P = pu(Y — 1) (u is the internal 
energy per unit mass), and the fluid is initially given the azimuthal velocity required for 
centrifugal equilibrium. The self-gravity of the disk is neglected, and the black hole produces 
a Newtonian point-mass potential $bh — —GM^-n/r. Accretion is modeled by an absorbing 
boundary at rs c h — 2GM^yi/c 2 . In practice, only the inner regions of the 3D disk are mapped 
onto two dimensions to start the calculation (with < \z\ < 200 km and < r < 400 km). 

The above simplifications clearly do not correspond to the physical conditions one would 
encounter in such a disk, but we have chosen to use them as a first approximation to the 
solution for several reasons. In the first place, the 3D dynamical calculations that gave rise 
to the disks we now use as initial conditions were fully Newtonian (except for the inclusion of 
gravitational radiation reaction terms), particularly regarding the form of the gravitational 
potential produced by the black hole. For consistency we have thus kept this form. In the 
future we will explore situations in which the black hole potential is given by a form which 
reproduces the effects of a marginally stable orbit (e.g. that of Paczyhski & Wiita (1980)). 
Second, the equation of state was likewise very simple in the initial calculations, and we have 
maintained it in the present work. The only change we apply is that at t — the value of the 
index T has been lowered for some runs (see Table 1) from 2 to 4/3. We will also improve 
on the form of the equation of state in future simulations. Third, the final conditions of the 
3D calculations show that the self-gravity of the disk is not too important in determining 
its structure (it is pressure support that makes the rotation curve deviate from what would 
be expected for point masses orbiting a massive body, see e.g. Figure 8b in L01). 

The second important new ingredient (the first one being the assumption of azimuthal 
symmetry) in the present simulations is the transport of angular momentum by viscosity, 
which is modeled by including all the terms derived from the stress tensor t a p in the equations 
of motion (we use the formalism developed by Flebbe et al. (1994)). We have: 



where the commas indicate partial derivatives. The last term arises from the presence of 
artificial viscosity (see below). The tensor responsible for the viscous force is t a p, given by: 




(1) 



t a p = r){v a ,f3 + Vp t , 



■ a 




(2) 
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where rj is the dynamical viscosity coefficient, and the shear is: 



2 

0"a/3 = Va,P + Vp,a ~ (3) 

The energy dissipation per unit mass is: 

r^ds 7] ( 'ds\ 

To perform a dynamical calculation, at this point one needs to fix a prescription for 
the viscosity. We have chosen to parametrize its magnitude by using an a-law, setting 
rj = apc 2 /Vt k , where c s = (TP/p) 1 ^ 2 is the adiabatic sound speed and Q k is the Keplerian 
angular velocity. This formulation of the viscosity has been used before in time-dependent 
studies of accretion flows in two and three dimensions (Igumenshchev, Chen & Abramowicz 
1996; Igumenshchev & Abramowicz 1999, 2000; Igumenshchev, Abramowicz & Narayan 2000; 
McKinney & Gammie 2002). We note that the simulations reported here are followed for a 
timescale that is much longer than the viscous one, and thus solutions that are qualitatively 
steady are obtained in all cases (they cannot be steady quantitatively simply because there 
is no external agent feeding the disk with matter as the simulation progresses). 

The implementation of SPH for the present 2D azimuthally symmetric calculations 
requires the writing of the above expressions in cylindrical coordinates. For completeness, 
we give here the full set of equations for momentum and energy: 

dv r ldP GM BU r 1 (dt„ t„ dt rz \ ( dv r \ 

+ -A — + — + — ) + ( — ) ( 5 ) 

art 



dt p Or R 3 p \ Or r dz J \ dt 



1 / 2t r dt^ r dt^z 



dt p \ r dr dz ) ^ 

dv z ldP GM BU z 1 fdt zr t zr dt zz \ fdv z \ 

dt ~ pdz R 3 p\dr r dz J \ dt ) art 1 j 

du ( P\ „ . ds , , 

-JV-u + T— (8) 



dt VP/ dt 
Here R = \Jr 2 + z 2 is the distance to the origin, and the artificial viscosity terms are com- 
puted using the prescription of Balsara (1995) to minimize the effects of artificial shear on 
the evolution of the disk. The explicit forms of the acceleration and the energy dissipation 
for a given SPH particle are: 

(§) =-E^iw% (9) 

\ / i,art j-^i 
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and 



V al J i,art Z j^i 

where W is the SPH interpolation kernel given by Monaghan & Lattanzio (1985), adapted 
to azimuthal symmetry: 



10 1 



1 - 3 {r/hf /2 + 3 (r/hf /4, < r/h < 1, 



^•^=7^(2^)1 P-^M l^r/k<2, (11) 



We have 



% = (5 + 5) (-«^« + ^4), (12) 



where 



{Vi-VjHVi-rj) h+fj 



(Vi - ■ (r, - rA < 



Mm-r-,| 2 /^)+r,2 2ci . , yvi >jj 

0, (Vi - Vj) ■ (Vi - Tj) > 



and fi is the form-function for particle i defined by 

fi = 





V • v . 


i 


V • v\ 


i + \ 


V x d 





The factor 77' ~ 10~ 4 in the denominator prevents numerical divergences. The smoothing 
length and sound speed for particle i are denoted by hi and q respectively, and and fib 
are constants of order unity (we use a& = /3 6 = T/2). The number of neighbors for the 
hydrodynamical interpolation is kept fixed at v — 24. 

Artificial viscosity is used in SPH to model the presence of shocks and avoid particle 
inter-penetration. The forms initially employed have been substantially improved over the 
years, with an accompanying reduction in spurious numerical effects. Most notably, the one 
we use here is such that the shear viscosity is suppressed when the compression in the fluid is 
low and the vorticity is high (| V xv\ 3> | V-#|, as is the case in accretion disks with differential 
rotation), but remains in effect if compression dominates in the flow (|V • v\ ^> |V x v\). 

We are able to monitor just how important the dissipation due to artificial viscosity is, 
compared with that arising from the stress tensor t a p (see equation 4). We employed two 
different values of the viscosity parameter a for the runs shown here (see Table l):a = 0.1 
and a = 0.01. For both cases, the energy dissipated through artificial viscosity is most 
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important at early times (t < 10 ms), and quickly decreases relative to the terms generated 
through t a p (see below, §3). 

We have tested our code to ensure that the viscous terms are computed correctly. The 
most simple test is to calculate the spreading of a ring of gas due to viscosity (with no vertical 
structure) orbiting a point mass, when pressure effects and self-gravity are neglected. This 
problem has an analytic solution in terms of Bessel functions (Pringle 1981), which can be 
easily compared with a numerical solution. Our code reproduces this calculation accurately, 
with the ring spreading radially on the correct timescale. 

We have also performed a second class of tests, allowing for vertical structure in az- 
imuthal symmetry, studying the evolution of thick tori, where pressure effects are important 
and self-gravity is neglected. Two different types of calculations were performed. In the 
first, a torus in hydrostatic equilibrium around a black hole, with an initially specified dis- 
tribution of specific angular momentum, was allowed to evolve, switching the viscosity on 
at the start of the calculation (Igumenshchev, Chen & Abramowicz 1996). In the second, 
a continuous flow of matter was injected at a large radius, and a disk was formed on the 
viscous timescale (Igumenshchev & Abramowicz 1999). For each of these, several values of 
the a parameter were used, ranging from 1CT 3 to 10 _1 . We found excellent agreement with 
previously published results, both qualitatively and quantitatively (Lee 2001b). 

Finally, we assume that all the energy dissipated in the flow is radiated away in neutrinos, 
and thus in this respect, the flow is an NDAF (Neutrino Dominated Accretion Flow), as done 
by Narayan, Piran & Kumar (2001) and Kohri & Mineshige (2002). This is accomplished 
by directly removing from the fluid the internal energy associated with the first term in 
equation (4). The neutrino luminosity is then given by L v = J ri<j a p<7 a p/2p dm. The energy 
dissipated by artificial viscosity remains in the disk as internal energy. This is another aspect 
that clearly needs to be treated in a more realistic manner and will be dealt with in future 
work. 



3. Results 

3.1. Structure and fluid circulation in the disks 

The choice of parameters shown in Table 1 for the simulations allows us to explore 
the effect of numerical resolution on the results (runs D and E, and runs F and G differ 
only in the initial number of particles), as well as that of the magnitude of the viscosity, the 
compressibility of the gas in the disk, and the initial mass of the black hole. For a lower black 
hole mass Mbh the accretion disk is initially more massive, and the spin angular momentum 
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of the black hole is larger. This is a reflection of the previous binary evolution and tidal 
disruption that led to the formation of the accretion disk. 

The initial conditions we have used for the accretion disks clearly do not correspond 
to a simple equilibrium configuration, since they are derived from three-dimensional dy- 
namical calculations. In our case, the disks are initially hot (the internal energy is about 
10 MeV/nucleon, see Figure 8 in L01), and thus pressure support is important, making them 
quite thick (H/R w 0.5). This energy reservoir is not entirely lost, since only the energy dis- 
sipated through the physical viscosity is radiated away in neutrinos. Thus the disks remain 
thick throughout the dynamical evolution. 

All the dynamical calculations are qualitatively similar. We show snapshots of the 
density and velocity field in the disk at various times for runs E (a = 0.1) and G (a = 
0.01) in Figures 1 and 2. At the start of the calculation, the disk flattens slightly, on a 
vertical free-fall timescale tff ~ ^t 3 /GMbr m 1 — 2 ms, since it is not in strict hydrostatic 
equilibrium at t — 0. The initial panels (at t — 10 ms) show the formation of a weak shock 
at r ~ 140 km, formed when the outer regions of the torus move radially inwards. This 
front moves subsequently outwards, and has practically left the inner region of the disk by 
t = 20 ms. Thereafter accretion proceeds onto the black hole in a characteristic pattern, 
which is mainly dependent on the value of a. 

The simulations use a variable smoothing length h to ensure an accurate interpolation 
of the fluid variables, and so the spatial resolution (which is essentially equal to h), is also 
variable. It is smallest in the regions of highest density (in the midplane of the disk). At 
t — 0.1 s, the smallest scales that can be resolved are of order 2 x 10 4 cm=200 m at r ~ 30 km, 
and 3 x 10 5 cm=3 km close to the surface of the disk (where p ~ 10 9 g cm" 3 for the high- 
resolution runs E and G; in two dimensions the smoothing length scales with the number of 
particles as h oc iV -1 / 2 ). The circulation patterns seen in Figure 2 are much larger than the 
spatial resolution and thus are clearly resolved. The runs with a lower number of particles 
have a slightly poorer resolution, but nevertheless resolve the same global features clearly in 
terms of accretion rate, spatial extension of the disk and circulation patterns. 

For a high viscosity, a large-scale circulation is established, with z ~ r. There is 
outflow in the equatorial plane at large distances from the black hole, and accretion takes 
place only along the surface of the disk and in the equatorial plane quite close to the black 
hole (r ~ 50 km). The radial velocities are subsonic, and the turnaround in v r occurs in 
large eddies close to the surface of the disk, when the gas is moving away from the equatorial 
plane. The disk slowly spreads in the radial direction, with an accompanying drop in density 
and accretion rate (at t = 0.2 s, M has dropped by about one order of magnitude compared 
to its value at 10 ms, see also Figure 7a below). 
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For the low viscosity case, the flow is more complicated. Accretion onto the black hole 
still occurs only along the surface of the disk and in the equatorial plane at small values of r. 
However, the circulation patterns are present on smaller scales than before, with z < r. The 
flow is qualitatively steady after approximately 20 ms, but the eddies persist throughout the 
simulation. The radial spread of the disk is less pronounced than before, due to the lower 
value of a, and thus higher densities are maintained for a longer period. The accretion rate 
onto the black hole is about one order of magnitude lower (for t < 20 ms) than for the high 
viscosity runs, and decreases more slowly. 

The circulation patterns described above can be clearly related to the distribution of 
specific angular momentum / in the disk. We show in Figure 3 the contours of / for runs E 
and G at t = 40 ms. As one moves from the equator (where v r > 0) to the surface of the 
disk at a given value of r, / decreases until the fluid turns around and flows inward. 

For both cases, the magnitude of the velocities (radial and vertical) in the disk decreases 
in time, with the circulation patterns becoming less intense in that respect. Figure 4 shows 
views of the disks for runs E and G at t = 0.11 s on a larger scale. 

The mass accretion rate onto the black hole is clearly sensitive to the value of a, and to 
a lesser extent, on the adiabatic index T, as can be seen in Figure 7a. Runs E and G differ 
only in the value of a, and clearly lowering the viscosity leads to less vigorous transport of 
angular momentum, and thus to a reduced M B h- Runs E and C differ only in the value of 
T. The less compressible case (run C, with T = 2) has a broader spatial mass distribution 
and reaches lower peak densities. Thus slightly more mass is initially closer to the accretor 
and easily absorbed, giving a larger Mbh, and subsequently a faster drop (by t = 0.2 s, the 
disk mass is lower in run C than in run E, see Table 2). Finally, run B differs from run E in 
the initial binary system that gave rise to the accretion disk, with a larger Mbh and smaller 
Mdisk, which is responsible for the lower accretion rate. 

Meridional circulation patterns like the ones seen in the present simulations have been 
studied before, both analytically and numerically. Solving for the vertical structure of an 
accretion disk and using an a viscosity, Urpin (1984) found that the radial flow can change 
direction in the midplane of the disk, with v r > 0. Kita (1995) and Kluzniak & Kita (2000) 
assumed a steady state in a polytropic a-disk and found solutions which exhibited outflow 
along the equator, with inflow along the disk surface and in the equator close to the accretor, 
as seen here. In their study it was the value of a that fixed the equatorial distance at which 
the flow turned from inflow to outflow (the stagnation radius). In both cases, the authors 
assumed that the disk cools efficiently, with the dissipated energy being radiated away. 

To our knowledge, numerical work has been performed so far mainly with Eulerian codes, 
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in which boundary conditions need to be specified. Kley & Lin (1992) obtained outflow in 
the midplane of the disk, and inflow along the surface, as in the solution found by Urpin 
(1984), but the use of height-averaged boundary conditions enforced artificial circulation 
at the edges of the computational domain, thus altering the flow structure substantially. In 
other cases, a free outflow condition is used, meaning that gas which leaves the computational 
domain cannot re-enter it at a later time. These calculations — which, contrary to what 
we have assumed, suppose that the dissipated energy is advected with the flow — have also 
revealed circulation patterns and inflows/outflows with morphologies that depend on the 
magnitude of the viscosity (Igumenshchev & Abramowicz 1999; Stone, Pringle & Begelman 
1999; McKinney & Gammie 2002). 

The use of SPH is clearly an advantage in this sense, since there are no boundary 
conditions at all and the modeling of the fluid is not restricted to a particular region in 
space. 

3.2. Dynamical stability 

There are several instabilities that can affect massive accretion disks dynamically and 
shorten their lifetimes considerably, compared with the viscous timescale. One of these is 
the "runaway radial" instability, first discussed in the context of massive accretion disks at 
the centers of galaxies (Abramowicz et al. 1983). A number of different effects can work 
to enhance or suppress this phenomenon, namely (i) the spin of the black hole, (ii) the 
distribution of angular momentum in the disk, given by / oc r p , (iii) the effects of general 
relativity and (iv) the self-gravity of the disk. The first two (large Kerr parameter and a 
high value of p) tend to suppress the instability (Wilson 1984; Daigne & Mochkovitch 1997; 
Abramowicz, Karas & Lanza 1998; Masuda, Nishida & Eriguchi 1998; Lu et al. 2000), while 
the others tend to enhance it (Nishida et al. 1996; Nishida & Eriguchi 1996; Masuda, Nishida 
& Eriguchi 1998; Font & Daigne 2002). The numerical approach used here allows us only to 
comment on (ii). The distribution of angular momentum in the disk can be approximated 
by a power law in the three-dimensional SPH coalescence simulations, with / oc r - 45 , i.e. 
quite close to Keplerian (see Figures 10b and 8b in Lee (2000) and L01 respectively). This 
does not vary appreciably as the calculations progress, despite the circulation of the fluid in 
the disk described above in §3.1, and thus we find that the disk is stable in this respect. 

Additionally, disks can be unstable to axisymmetric perturbations depending on their 
surface density and temperature (measured through the sound speed c s ). This is known as 
Toomrc's stability criterion, and can be quantified through Q = c s k/itGY,, where k is the 
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local epicyclic frequency and X is the mass surface density: 

/oo 
p(r,z)dz. (13) 
-oo 

For Q > 1 the disk is stable to axisymmetric perturbations, and unstable otherwise. This 
condition tells us that a dense, cold disk is more likely to be unstable than a hot, thick 
one (such as the ones treated here). We find that for all our runs, Q > 1. The profile 
of Q is such that at a given instant in time, there is a minimum Q m in at a certain radius 
(close to the maximum in density, see Figure 5). The lowest value of Q m in occurs at t — 0, 
where (rq, Q m i n ) ~ (50 km, 5). As the disk evolves, the density drops and both parameters 
increase. By t = 0.2 s, we typically find (vq, Q m in) ~ (150 km, 10). 



3.3. Relative importance of artificial viscosity 

The artificial viscosity mentioned above in §2 ideally should not affect the accretion 
flow substantially, otherwise the results would be difficult to assess. As a measure of the 
effect the terms arising from it have on the disk, we plot in Figure 6 the dissipated energy 
as a function of time, separated into the components arising from t a p and those coming 
from artificial viscosity, for runs E and G (the former is in fact what we have termed the 
neutrino luminosity L u above in §2). It is at early times (t < 10 ms) that artificial viscosity 
is important. This is not surprising, since that is when the initial transient appears in the 
disk and the shock forms at r ~ 140 km. Thus the relative importance of the artificial 
viscosity merely means that it is working as one would expect. After this transient dies 
out, the dissipation rate becomes quickly dominated by the physical viscosity. It is also 
apparent that only for low values of a (run G with a = 0.01) are the two terms comparable 
in magnitude in the initial stages. We are thus confident that our results are not affected 
significantly by the presence of artificial viscosity. 



4. Summary and discussion 

We have computed the dynamical evolution of massive accretion disks around stellar- 
mass black holes in two dimensions (with azimuthal symmetry), formed through the tidal dis- 
ruption of a neutron star by a black hole in a close binary (M B h ~ 4M and M disk « O.3M ). 
Our initial conditions are taken from the final state of three dimensional hydrodynamical 
calculations of the coalescence process, by averaging in the azimuthal direction. We use 
Newtonian physics, an ideal gas equation of state, and solve the equations of viscous hy- 
drodynamics assuming an a law for the viscosity coefficient. All the energy dissipated by 
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the physical viscosity is radiated away (in neutrinos). The time evolution is followed for 
0.2 seconds. 

We find that after an initial transient of numerical origin, stemming from the fact 
that the 3D torus does not exhibit strict azimuthal symmetry due to the highly dynamical 
merging process (see Figures 2 and 7 in L01), the disk settles to a qualitatively steady state. 
Meridional circulations are promptly established, whose structure depends mainly on the 
value of the a parameter. Most strikingly, there is an important motion of fluid from the 
inner regions of the disks to large radii, along the equatorial plane. The flow is directed 
towards the accreting black hole along the surface of the disk and in the equatorial region 
at small radii. The disks remain thick {H/R ~ 0.5) throughout the dynamical evolution, 
due to their large internal energy, with accretion rates on the order of one solar mass per 
second. The maximum densities decrease during our calculations, as there is no external 
agent feeding the disks, but remain at ~ 10 12 g cm~ 3 , with corresponding internal energy 
densities ~ few x 10 30 erg cm" 3 . 

We stress that the evolution of accretion disks such as these should be studied with time- 
dependent models, since the system is clearly not in a steady state, even from its inception. 
The circulation pattern seen in Figure 2 would certainly lengthen their lifetimes by moving 
matter to larger radii continuously, an effect that would otherwise be omitted. To illustrate 
how important this can be, we considered the total radial mass flux in the disk, composed 
of two parts at any given value of the radial coordinate r, M(r) = M in (r) + M out (r) where 

M in = 2nr / pv r dz, M out = 2nr / pv r dz, (14) 

Jv r <0 J v r >0 

restricted to regions in which v r < for M in and v r > for M out respectively. By definition, 
M in < and M out > 0. If there were no circulation in the disk and all matter moved 
radially inward one would have M out = and M = Mi n < 0. A measure of how important 
the circulations are, and how much they would lengthen the lifetime of the disks can be 
obtained by calculating the fraction of the gas flowing radially in the disk that is actually 
moving toward the accretor, i.e. |M i7l |/(|M in | + |M out |). In the inner regions of the disks 
this ratio tends to unity, as can be seen from Figure 2. It decreases rapidly at larger radii, 
reaching about 1/3 for run E and 1/10 for run G at r ~ 100 km midway through the 
simulations (at t = 0.1 s). Thus most of the radial flow of gas actually cancels out in the 
circulations, with a residual amount left over moving toward the black hole. 

We now turn to the implications our calculations might have on models for the pro- 
duction of cosmological gamma ray bursts from coalescing compact binaries. The two main 
forms of energy release from the disk we consider are i) neutrino emission and ii) MHD flow, 
either through the Blandford-Znajek effect (Blandford & Znajek 1977) or by means of a 
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magnetized wind. 

In the first case, we make a rough estimate for L v and E v from the dissipated energy 
because of viscosity (see Table 2), as mentioned above. We have made an estimate only for 
the total neutrino luminosity L u , and not for the annihilation luminosity L uV , which would 
determine if a relativistic fireball could be launched or not. The calculation of L vV requires 
the use of a more realistic equation of state, which we will explore in future work. For the 
time being, we note that, regardless of the efficiency of energy conversion from neutrino 
luminosity to annihilation luminosity — which could be quite low, on the order of 1 per 
cent or less (Ruffert & Janka 1999; Popham, Woosley & Fryer 1999) — it appears that the 
time-dependence of L„, which follows that of M B h (see Figure 7a and Table 2), is such 
that neutrinos could only be responsible for a very short, almost impulsive energy release 
(L v oc t~ 5//4 , so L uV oc L 2 oc t~ 5 / 2 ), and thus would be unable to power a burst lasting several 
tenths of a second or more. Of course this does not mean that it would have a negligible 
impact on the structure of the burst itself. We note furthermore that in the detailed 3D 
calculations done by Ruffert & Janka (1999) for the evolution of thick disks following the 
coalescence of two neutron stars, the "neutrinosphere" is quite close to an isodensity surface 
at p = 10 11 g cm~ 3 , which is lower than the maximum densities present in the disks computed 
here, even at late times (see Figures 1 and 4). Rosswog & Davies (2002) have also performed 
3D calculations of binary neutron star coalescence taking into account neutrino emission 
and scattering processes, finding as well that in the regions of highest density the material 
is opaque, in their case mainly because of scattering off heavy nuclei. Thus it is clear that a 
complete picture must include an appropriate formulation of neutrino transport. 

For the magnetic-dominated case we must make some assumptions, since our simula- 
tions do not incorporate the effects of magnetic fields explicitly. For the field to be able 
to extract the binding energy of the torus, it should be anchored to it, and we assume its 
magnitude is directly related to the internal energy in the gas. For this purpose we show 
in Table 2 the internal energy density pc 2 at t — 0.1 s in the inner regions of the disk (at 
r = 1.25rsch (~ 15 km for runs A and B, and m 20 km for runs C through G) in the equa- 
torial plane. It is of order 10 30 erg cm -3 , and even larger for the case with low viscosity 
(a = 0.01). From this we compute an estimate for the Blandford-Znajek luminosity as 

L ^^(M(^oy^ s " (i5) 

where a ~ 0.3 is the Kerr parameter of the black hole and the magnitude of the magnetic 
field is computed using B 2 /8n = pc 2 s (this gives B 10 16 G in all cases). This is clearly 
the most optimistic scenario concerning energy release, in that it assumes that the magnetic 
field strength is at the equipartition value. The time evolution of L BZ is shown in Figure 7b 
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for several runs. For a larger viscosity, the gas in the disk drains into the black hole on a 
shorter timescale, and thus the drop in L BZ is much faster than for a low value of a (in run 
G, L BZ 5 x 10 51 erg s -1 is practically constant). 

In principle, energy can be extracted over many dynamical timescales if magnetic fields 
can tap the rotational energy of the accretion disk and of the black hole. Here L BZ oc (t/t ) q ' 
is the intrinsic luminosity of the central engine measured in the fixed frame, where —5/4 < 
q' < (see Figure 7b). For such a continuously fed fireball, a forward shock propagating 
into the external medium and a reverse shock moving back into the relativistic outflow 
will coexist on either side of the contact discontinuity. The latter may persist as long as 
a significant level of energy injection is maintained. The differential energy conservation 
relation for the self-similar blast wave can be written as dE/dt = L BZ (t/t ) q ' — n'(E/t) 
(Cohen, Piran & Sari 1998), where the first term denotes steady energy injection, and 
the second accounts for possible radiative energy losses of the fireball, with q' and k' being 
dimensionless constants. In the adiabatic case (no radiative losses) k' — and the Lorentz 
factor evolves as V 2 oc t~ 3 (Blandford & McKee 1976). For q' ^ — 1 — k', an analytical 
solution can be found: E = [L BZ /(q' + «' + l)](t/t ) 9 + E imp (t/t )~ K for t > t , where 
E imp = E vV describes the impulsive energy input. Here to is the characteristic timescale for 
the formation of a self-similar solution, which is roughly equal to the time for the external 
shock to start to decelerate. For t > t , the bulk Lorentz factor of the fireball scales with 
time as T 2 oc t~ m , with m and k' connected by k' = m — 3, and m = 3 for an adiabatic 
blast wave expanding in a constant density medium. In the observer frame, the arrival time 
at the detector T is related to that in the fixed (laboratory) frame t by dT — (1 — (3)dt, 
and T = J Q t (2T 2 )- 1 dt ^ t/[2(m + l)r 2 ] (Fenimore et al. 1999). The differential energy 
conservation relation in the observer frame is now given by dE/dT = C B z{T/To) q — n(E/t), 
and the integrated relation is 

B = ^TTS) ,;r + ^S)""- r>T °' (16 » 

Here C BZ = 2T 2 (t )L BZ , and q = (q' — m)/(m+ 1), k = k' / (m + 1). Since k + q + 1 = 
(q' + k' + l)/(m + 1), the comparisons between q and — k — 1 in the following discussion 
are equivalent to the comparisons between q' and —n' — 1. At different times, the total 
energy of the blast wave may be dominated either by a continuous injection term or by the 
initial impulsive term (Zhang & Meszaros 2001, and references therein). Which of the two is 
dominant at a particular observation time T depends both on the values of the two indices 
(q and — k — 1), and on the relative values of C BZ and E vV . If q < —1 — k, the impulsive term 
always dominates since the first term is negative. For an adiabatic blast wave expanding 
in a constant density medium (i.e. m = 3 and k — k' — 0), this condition yields q' < — 1. 
If q > — 1 — k, the continuous term will eventually dominate (i.e. E vV C BZ T ) over the 
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second term after a critical time 

T c = T x max 



1, (/c + g + 1) 



CbzTq 



(17) 



This is the case considered in many pulsar central engines (Usov 1992, 1994; Dai & Lu 1998). 
We then have T c T and the fireball is completely analogous to the impulsive regime with 
-E'total ~ CbzTq. 



If initially E vV ^> CbzTq, the critical T c after which the continuous injection becomes 
dominant could be much longer than T , exerting a noticeable influence on the GRB after- 
glow. We expect this to be the case either when the magnetic energy is below equipartition 
(pCg > B 2 /8n) or when the timescale for the formation of a self-similar solution is small 
(T <C T c ). Because of the strong pinch that develops, a narrow jet that delivers its thrust 
in a narrow solid angle, Qbz, m ay be a common ingredient of strong rotating magnetic 
fields (see Meier, Koide & Uchida (2001) for a recent review). We thus generally expect 
the magnetic outflow to be much more collimated than that produced by vv annihilation 
(Qbz ^ Qvv) and the magnetic luminosity to be the dominant contribution. Even in this 
case, the impulsive term E vV may be responsible for either creating a cavity before the mag- 
netized wind expands or for precursor emission. The detection of, or strong upper limits 
on, such features would provide constraints on the burst progenitor and on magnetar-like 
central engine models. 

An external shock can occur at much larger radii and over a much longer timescale than 
in usual afterglows, if the environment has a very low density. This may be the case for 
GRBs arising from compact binary mergers that are ejected from the host galaxy into an 
external medium that is much less dense than the ISM assumed for usual models (where the 
particle density is n ~ 0.1 — 1 cm -3 ). It is commonly assumed that compact mergers occur 
outside the host galaxies because of the long inspiral, due to the emission of gravitational 
waves. We note that recent calculations of evolutionary tracks of binary systems containing 
massive stars show that the merger events can occur on much shorter timescales and thus 
still within the host galaxy, because the initial binary separation is small (Belczynski, Bulik 
& Kalogera 2002). 

Our simulations would indicate that the central engine survives the initial, violent event 
that created it, and that it possesses enough energy to account for the energetics and du- 
rations of short GRBs. However, they clearly cannot tackle directly other relevant issues, 
mainly related to the evolution of the magnetic field, and its influence on the dynamics. 
Magnetic instabilities could make the disk lifetime much shorter by effectively increasing 
the viscosity. The amplification of the magnetic field may be self-limiting due to magnetic 
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stress, which would cause disk flaring. The properties of the expected variability depend 
strongly on the details of the configuration of the disk corona generated by the magnetic 
field, which is removed from the disk by flux buoyancy (Narayan, Paczyhski & Piran 1992; 
Tout & Pringle 1996; Kluzniak & Ruderman 1998). 

It is a pleasure to acknowledge many helpful conversations with W. Kluzniak, D. Laz- 
zati, G. Ogilvie, D. Price, M. J. Rees, S. Rosswog and V. Usov. We thank the referee for 
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hospitality. Financial support for this work was provided in part by CONACyT (27987E), 
PAPIIT (IN- 11 0600), SEP and the Royal Society. 
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Fig. 1. — Logarithmic density contours (equally spaced and labeled in g cm~ 3 ) at t — 10 ms 
(top left), t — 20 ms (top right), t — 30 ms (bottom left) and t — 40 ms (bottom right) for 
runs E and G. In each panel, the left half (r < 0) corresponds to run E (a = 0.1) while the 
right half (r > 0) corresponds to run G (a = 0.01). 
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Fig. 2. — Same as Figure 1 but showing the velocity field. The equatorial outflow at large 
radii, and the inflow along the surface of the disk can be clearly seen, particularly for run E 
(left half of each panel). Note the small-scale eddies present in run G in every panel. The 
strength of the circulation diminishes in time, but does not die out. The velocity field is only 
shown if p > 6 x 10 8 g cm -3 for clarity (in these panels the edge of the disk, as given by the 
SPH particle distribution extends to lower densities, by about one order of magnitude, see 
Figure 1). The vector in the top left corner of each panel has magnitude v — 5 x 10 8 cm s _1 . 
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Fig. 3. — Contours of specific angular momentum / for runs E (a = 0.1) and G (a = 0.01) 
at t = 40 ms (corresponding to the last panel in Figures 1 and 2). The left half of the panel 
(r < 0) corresponds to run E, and the right half (r > 0) to run G. The contours are equally 
spaced and the labels correspond to //10 16 cm 2 s -1 . 
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Fig. 4. — Logarithmic density contours (equally spaced and labeled in g cm -3 ), and velocity 
field for run E (a) and G (b) at t — 0.11 s. Note the different scales on the axes for each 
panel. The longest arrows correspond to v — 2.9x 8 cm s -1 for (a) and v = 10 8 cm s" 1 for 
(b). 
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Fig. 5. — Toomre stability parameter Q as a function of the radial coordinate r for runs E 
(solid line, a = 0.1) and G (dashed line, a = 0.01), at t — 0.1 s. 
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Fig. 6. — Energy dissipation rate as a function of time for runs E and G arising from the 
stress tensor t a p (solid lines, L u ), and from artificial viscosity (dashed lines). 
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Fig. 7. — (a) Accretion rate (in solar masses per second) onto the black hole for runs B, C, 
E and G. (b) Blandford-Znajek luminosity Lbz for the same runs as shown in (a). 
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Table 1. Initial conditions for the accretion disks. 



Run 




r 


a 


M BH /M b 


M disk /M Q h 


Arc 


A 


0.31 


2.0 


0.1 


5.57 


0.266 


21,609 


B 


0.31 


4/3 


0.1 


5.57 


0.266 


21,609 


C 


0.50 


2.0 


0.1 


3.85 


0.308 


19,772 


D 


0.50 


4/3 


0.1 


3.85 


0.308 


19,772 


E 


0.50 


4/3 


0.1 


3.85 


0.308 


79,489 


F 


0.50 


4/3 


0.01 


3.85 


0.308 


19,772 


G 


0.50 


4/3 


0.01 


3.85 


0.308 


39,658 



a Mass ratio g& = Mns/Mbh of the original binary system 
used in the three-dimensional simulation. 

b Values at the start of the two-dimensional calculation of 
the accretion disk evolution. 

c Number of SPH particles at the start of the two- 
dimensional calculation. 
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Table 2. Accretion disk parameters during the dynamical evolution. 



Run 


M/(M Q /s) a 


pc 2 / (erg cm 3 ) a 


L^crgs- 1 )- 


£„(erg) b 


L BZ (crgs x ) a 


Ebz (erg) b 


M dl3k /M Q h 


A 


0.33 


1.8-10 30 


3.25-10 52 


1.38-10 52 


1.00-10 51 


2.94-10 50 


0.128 


B 


0.31 


2.1-10 30 


3.61-10 52 


1.08-10 52 


1.75-10 51 


4.05-10 50 


0.151 


C 


0.36 


3.5-10 30 


4.00-10 52 


1.79-10 52 


2.00- 10 51 


6.47-10 50 


0.137 


D 


0.48 


5.0-10 30 


5.00-10 52 


1.52-10 52 


2.50- 10 51 


6.25-10 50 


0.158 


E 


0.47 


5.3-10 30 


4.65-10 52 


1.57-10 52 


2.95- 10 51 


7.41-10 50 


0.160 


F 


0.13 


12.5-10 30 


1.86-10 52 


0.44-10 52 


6.00-10 51 


1.17-10 51 


0.254 


G 


0.15 


12.0-10 30 


1.84-10 52 


0.44-10 52 


5.65-10 51 


1.15-10 51 


0.255 



a Valucs arc given at t = 0.1 s. 
b Values are given at t = 0.2 s. 



